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Electrons in certain two-dimensional crystals possess a pseudospin degree of freedom associated 
with the existence of two inequivalent valleys in the Brillouin zone. If, as in monolayer M 0 S 2 , 
inversion symmetry is broken and time-reversal symmetry is present, equal and opposite amounts of 
k- space Berry curvature accumulate in each of the two valleys. This is conveniently quantified by the 
integral of the Berry curvature over a single valley - the valley Hall conductivity. We generalize this 
definition to include contributions from disorder described with the supercell approach, by mapping 
(” unfolding”) the Berry curvature from the folded Brillouin zone of the disordered supercell onto 
the normal Brillouin zone of the pristine crystal, and then averaging over several realizations of 
disorder. We use this scheme to study from first-principles the effect of sulfur vacancies on the 
valley Hall conductivity of monolayer M 0 S 2 . In dirty samples the intrinsic valley Hall conductivity 
receives gating-dependent corrections that are only weakly dependent on the impurity concentration, 
consistent with side-jump scattering and the unfolded Berry curvature can be interpreted as a k- 
space resolved side-jump. At low impurity concentrations skew scattering dominates, leading to 
a divergent valley Hall conductivity in the clean limit. The implications for the recently-observed 
photoinduced anomalous Hall effect are discussed. 

PACS numbers: 71.15.Dx, 71.23.An, 72.10.Fk, 73.63.-b 


I. INTRODUCTION 

Monolayers of M0S2 and related transition-metal 
dichalcogenides (TMDs) have recently become the sub¬ 
ject of intense investigation, due in part to the possi¬ 
bility of manipulating the so-called “valley” degree of 
freedom. 1 These materials have the symmetry of a honey¬ 
comb structure with a staggered sublattice, thus lacking 
an inversion center. The bandstructure exhibits a direct 
gap at the two inequivalent valleys centered at the high- 
symmetry points K and K' — —K in the Brillouin zone 
(see Fig. 1), where the topmost valence bands are are 
primarily composed of transition-metal d states. 2 Time- 
reversal symmetry, which takes k into — k and therefore 
maps one valley onto the other, dictates that states in 
a given band at K and K' carry antiparallel angular 
momenta. This inspired Xiao et al. to propose using 
circularly-polarized light as a means of selectively excit¬ 
ing carriers from a particular valley. 3,4 The effect was 
rapidly confirmed experimentally, by demonstrating that 
excitation with circularly-polarized light results in polar¬ 
ized fluorescence. 5-7 

The broken inversion symmetry in monolayer M0S2 
induces a nonzero Berry curvature on the Bloch bands 
(in contrast, the Berry curvature vanishes identically for 
bilayer and bulk MoS 2 , both of which are centrosym- 
metric). The Berry curvature is defined in terms of the 
cell-periodic Bloch states as 

^ n,xy (k) — 21m {^k x ^n\^\&k y ^nk) 5 (1) 

n 

and it modifies the current response to an applied elec¬ 
tric field by adding an “anomalous velocity” term to the 


semiclassical equations of motion. 8 A well-known conse¬ 
quence is the anomalous Hall effect (AHE) in magnetic 
materials, where the Berry curvature is induced by bro¬ 
ken time-reversal symmetry. The intrinsic, part of the 
anomalous Hall conductivity (AHC) is given by the Bril¬ 
louin zone (BZ) integral of the Berry curvature summed 
over the occupied states, 8,9 

n*y(k) = Ti fnkftn,xy(k) (2) 

n 

a ° xv = ~~h J BZ ( 2 -)'' ’ (3) 

where f n k is the occupation factor and d is the dimen¬ 
sionality. Here, the superscript 0 denoted that it is the 
intrinsic part of the AHC. For d — 2 the AHC has units 
of conductance (S), and for d = 3 it has units of con¬ 
ductivity (S/cm). The broken inversion symmetry in 
monolayer MoS 2 induces a nonzero Berry curvature on 
the Bloch bands (in contrast, the Berry curvature van¬ 
ishes identically for bilayer and bulk MoS 2 , both of which 
are centrosymmetric). The Berry curvature is defined 
in terms of the cell-periodic Bloch states as Monolayer 
MoS 2 is nonmagnetic, and the presence of time-reversal 
symmetry implies the relation 8 

U^(-k) = -^(k). (4) 

Thus, equal and opposite amounts of Berry curvature 
accumulate in the two valleys, resulting in a cancellation 
of the valley Hall currents and a vanishing AHC. Time- 
reversal symmetry can however be broken by illuminating 
the sample with circularly-polarized light, leading to a 
photoinduced AHE. The valley-selective photoexcitation 
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FIG. 1. (Color online). Brillouin zone of monolayer M 0 S 2 , 
partitioned into two triangular valleys delimited by lines of 
vanishing Berry curvature, drawn in red. The two valleys are 
centered at the high-symmetry points K and K'. 


creates a carrier imbalance which in turn removes the ex¬ 
act cancellation between the Hall currents in the two val¬ 
leys. This so-called valley Hall effect was first discussed 
for graphene systems with broken inversion symmetry, 3 
and later for monolayer M 0 S 2. 4 The effect was subse¬ 
quently measured by Mak at al. in transistors of M 0 S 2 
monolayers. 10 

Compared to the conventional AHE in ferromagnetic 
metals, the theoretical modeling of the photoinduced 
AHE in TMDs poses the additional challenge that the 
AHC should in principle be calculated for a nonequi¬ 
librium photoexcited state, but to our knowledge, such 
a calculation has not yet been attempted. Instead, 
an approximate but more tractable approach is often 
used. 3,4,10 The idea is to introduce an auxiliary quan¬ 
tity crf yl which we will call the valley Hall conductivity 
(VHC), defined as the integral of the Berry curvature 
over a single valley domain in the BZ. For example, the 
intrinsic VHC of the valley centered at K in Fig. 1 is 

Vx V k (.£f) = f R H-n xy (k). ( 5 ) 

and similarly for the valley centered at K '. (The demar¬ 
cation of the two valley domains will be discussed further 
in Sec. II A. Equation (5) depends on the Fermi level ep 
through the occupation factors in Eq. (2).) The photoin¬ 
duced AHC a xy is then approximated by the sum of the 
VHCs in the two valleys, positing a Fermi-level shift Se 
between them to mimic the effect of the valley-selective 
photoexcitation, 

Sa xy {£ F , Se) = a^ y (e F + Se ) + cr^'(e F ) . (6) 

When Se = 0 the AHC vanishes, and a nonzero S<j xy 
appears when Se 0. This approach, also allows for a 
direct comparison with model calculations without con¬ 
sidering the details of how the carrier imbalance between 
valleys is generated. 

This will be the basic approach taken in the present 
work. We have suppressed the superscript 0 from this 
last equation to emphasize that it remains valid when 
the non-intrinsic contributions which we will now discuss 
are taken into account. 

Impurities are always present in real samples, and 
their extrinsic contributions to the photoinduced AHE in 


TMDs should be taken into account alongside the intrin¬ 
sic response described by Eqs. (5) and (6). This is well 
known in the context of the AHE in ferromagnetic met¬ 
als, where historically two types of extrinsic contributions 
have been considered - side jump and skew scattering . 9 In 
a simplified picture, the side-jump effect originates in the 
anomalous velocity that a wave packet may acquire as it 
moves through an impurity potential, while skew scatter¬ 
ing arises from the chiral part of a standard transition- 
rate expression. With some effort, both contributions 
can be incorporated into the semiclassical Boltzmann- 
transport framework. 11 

The correspondence between the semiclassical treat¬ 
ment of the AHC and a fully quantum-mechanical 
(Kubo-Streda) calculation based on a perturbative ex¬ 
pansion in powers of the disorder strength was carefully 
worked out in Ref. 12. It became clear from that anal¬ 
ysis that not all terms fall distinctly into either of the 
above physical interpretations of extrinsic contributions 
to the AHC, and for many purposes it is more practi¬ 
cal to base the distinction on the scaling with impu¬ 
rity concentration. 9 According to this viewpoint skew¬ 
scattering is defined as the part of the AHC which scales 
inversely with the impurity concentration, while the part 
which is independent of the impurity concentration has 
both intrinsic and side-jump components. Although the 
intrinsic contribution is sharply defined theoretically in 
terms of the electronic structure of the pristine crystal by 
Eq. (3), experimentally it is not known how to separate it 
from the side-jump part. Note that the anomalous Hall 
response of pristine samples at low temperatures is domi¬ 
nated by skew-scattering, with the intrinsic contribution 
only becoming significant in moderately resistive sam¬ 
ples (where it competes with side-jump scattering). This 
analysis, originally developed for the AHC in ferromag¬ 
netic metals, carries over to the VHC and photoinduced 
AHC in TMDs. 

It is well established that sulfur vacancies constitute 
the main source of disorder in M 0 S 2. 13 18 The forma¬ 
tion energies and thermodynamics of these defects have 
been thoroughly studied, 19 but their influence on trans¬ 
port and optical properties remains largely unexplored. 
Modeling the effects of disorder from first-principles is 
a challenging task in general, but there are notewor¬ 
thy examples where the AHC in ferromagnetic mate¬ 
rials has been calculated using the coherent potential 
approximation; 20 ' 25 also, an ab initio implementation of 
the side-jump contribution to the AHC has been car¬ 
ried out assuming scattering centers with delta-function 
potentials. 26 

In this work, we develop a computational scheme that 
allows us to include in a realistic manner the effect of 
impurities in the calculation of the VHC and of the pho¬ 
toinduced AHE in TMDs. In a first step, we perform 
several supercell calculations at the desired impurity con¬ 
centration, corresponding to different realizations of dis¬ 
order. In order to carry out the calculations efficiently 
while maintaining first-principles-like accuracy, we con- 
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struct effective Hamiltonians in a Wannier-function basis, 
starting from density functional theory calculations on 
smaller cells. 27 Recall that the definition of the VHC in 
Eq. (5) requires identifying the individual valley domains 
in the BZ where the Berry curvature is to be integrated. 
It is not clear a priori how to do so in the context of a 
supercell calculation, since the electronic states cannot 
be labeled by wavevectors in the normal BZ of Fig. 1. 
In order to overcome this difficulty, in a second step we 
use a “BZ unfolding” technique 28 to map the results of 
the super cell calculation onto the normal BZ of the pris¬ 
tine crystal. More precisely, we express the AHC of each 
disordered supercell configuration as an integral of the 
supercell Berry curvature over the folded BZ, and then 
unfold the Berry curvature onto the normal BZ accord¬ 
ing to the prescription of Ref. 29. Having done that, 
the VHC (including the contributions from disorder) can 
then be obtained by integrating the unfolded Berry curva¬ 
ture over a single valley domain in Fig. 1, and averaging 
the result over several realizations of disorder. 

We have used the above first-principles-based method¬ 
ology to study the influence of sulfur vacancies on the 
VHC of MoS 2 as well as the photoinduced AHC, which 
is the quantity measured in experiments. The calcu¬ 
lated VHC as a function of defect concentration was 
compared with model calculations where the valence and 
conduction-band edges in each valley are described by a 
massive Dirac Hamiltonian with a random distribution 
of delta-function scatterers. 12 

The paper is organized as follows. We start Sec. II by 
reviewing some basic features of the electronic structure 
of M 0 S 2 . We then evaluate the intrinsic VHC [Eq. (5)] 
and photoinduced AHC [Eq. (6)] for the massive Dirac 
Hamiltonian without disorder, and carry out the corre¬ 
sponding ab initio calculations for pristine M 0 S 2 . Our 
main results are presented in Sec. Ill, where disorder ef¬ 
fects are included in the calculation of the VHC, both for 
the massive Dirac Hamiltonian and for M 0 S 2 with sulfur 
vacancies. The two types of calculations are found to be 
in reasonable agreement, and we then proceed to calcu¬ 
late the photoinduced AHC for the disordered massive 
Dirac model as a function of gating voltage and Fermi- 
level shift, finding good agreement with the experimen¬ 
tal measurements. Our conclusions are summarized in 
Sec. IV, and the appendices present the details of the 
ab initio calculations, the BZ unfolding method, and the 
effective-Hamiltonian methodology. 


II. PRISTINE MOS 2 
A. Energy bands and Berry curvature 

Ab initio density-functional theory calculations were 
carried out for monolayer M 0 S 2 as described in Ap¬ 
pendix A 1. The calculated Kohn-Sham energy bands are 
shown in the upper panel of Fig. 2, color-coded by the 
spin component ( S z ) orthogonal to the layer. The mini¬ 


mum direct gap is situated at K and K' = —K, with a 
value of ~ 1.7 eV. Away from the time-reversal-invariant 
points T and M the spin degeneracy is split by the combi¬ 
nation of broken inversion symmetry and spin-orbit cou¬ 
pling (the degeneracy is actually protected along the en¬ 
tire T-M line by mirror symmetry). The two topmost 
valence bands exhibit a maximum spin-orbit splitting of 
~ 0.15 eV at K and IT', where S z is a good quantum 
number and Kramers-degenerate partners have opposite 
spin character: Sk^ = £k'1- 

The lower panel of Fig. 2 shows the Berry curvature 
summed over the valence bands, Eq. (2). In agreement 
with Eq. (4), Q xy is an odd function of k. Its magnitude 
peaks at the minimum-gap points K and iT', where the 
sign is dictated by the dominant contributions coming 
from the topmost valence bands. At K those bands are 
mainly composed of molybdenum d-states with mi — 2; 
according to the optical selection rules 7 those states can 
be excited with left-handed polarized light. The Berry 
curvature has a secondary peak between T and K\ there, 
a pair of lower-lying valence bands approaches the top¬ 
most ones, and also contributes significantly to the Berry 
curvature. 

It should be noted that because the Berry curvature is 
induced by the broken spatial inversion, it is not directly 
related to the spin-orbit splitting evident in Fig. 2; in 
fact, Q xy is practically unaltered if the spin-orbit inter¬ 
action is switched off. (This is in sharp contrast to the 
Berry curvature induced by broken time-reversal sym¬ 
metry in ferromagnetic metals, which vanishes in the ab¬ 
sence of spin-orbit coupling. 8,9 ) Likewise, the extrinsic 
scattering contributions do not rely on spin-orbit cou¬ 
pling. Thus one can use a spinless model such as the 
massive Dirac Hamiltonian of Sec. IIB 1 to describe the 
band edges and valley Berry curvature in M 0 S 2 ; spin is 
then accounted for by inserting a factor of two in the 
calculated Q xy . 

The Hamiltonian of monolayer M 0 S 2 is invariant under 
reflection across the vertical planes containing the lines 
that connect a Mo atom to the neighboring S atoms. The 
corresponding symmetry elements in reciprocal space are 
the T-M mirror lines. The Berry curvature transforms 
like a magnetic field in reciprocal space. 8 In particular, 
the component Vt xy = Dt z is odd under reflection across 
the T-M lines, and hence it vanishes along those lines, 
which form the boundaries between the two valleys in 
Fig. 1. This allows us to uniquely define the intrinsic 
VHC according to Eq. (5). 


B. Intrinsic Valley Hall conductivity 

1. Massive Dirac model 

The valence and conduction-band edges of a single val¬ 
ley of M 0 S 2 and related materials are often modeled by 
the massive Dirac Hamiltonian. 1,3,4 The Hamiltonian for 
the K valley in Fig. 2 reads 
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FIG. 2. (Color online). Top: Calculated bandstructure of 
monolayer M 0 S 2 . Energies are measured from the Fermi level, 
and the bands are color-coded by the spin expectation value 
(S z ), with red corresponding to spin up and blue to spin 
down. Bottom: Berry curvature summed over the valence 
bands [Eq. (2)]. The dashed line is the Berry curvature eval¬ 
uated from the two topmost valence bands. 


H( k) = hv(—k x a x + ky(jy ) + A o z , (7) 


where k = (k x , k y ) is measured from the valley center iF, 
(Ji are the Pauli matrices, and A is the mass parameter. 
The energy eigenvalues and Berry curvature are 


£±(k) = =b\/ A 2 + v 2 h 2 k 2 (8) 

A v 2 h 2 

C±,a:j/(k) = ±2( A 2 +v 2 fi 2 fc 2)3/2 ’ ^ 

where k 2 = k 2 -\-k 2 . In the case of the K' valley Eq. (8) 
remains unchanged, while Eq. (9) flips sign. 

If the Fermi level lies in the conduction band (£p > A), 
the VHC becomes 


a 


0 ,k _ 

xy 



2 h 2 h vA £p) 


f fVwl 

J\k\<k F Z7T J 


( 10 ) 


The valley and spin degrees of freedom can be included 
by coupling four copies of this model corresponding to 
spin and valley degrees of freedom. Referring to the 
band structure of M 0 S 2 we find that in the vicinity of 
the valleys the massive Dirac model provides a good fit 
if we take A = 0.86 eV. The velocity v is calculated as 
v 2 = A/m*, where m* = 0.4m e is the effective mass of 
the conduction band valley. 

For the ungated case where sp is in the gap, the model 
gives an intrinsic VHC of e 2 /2 h. The deviation from that 
result measures the contribution from lower lying valence 
bands and the crystal potential, which gives rise to non- 
hyperbolic bands away from the top of the valleys. It 


should be noted though, that a Chern insulating sys¬ 
tem with an inversion center will retain the exact value 
= e 2 / 2h when the crystal potential is included, be¬ 
cause symmetry implies fi(k) = fi(—k) and the topol¬ 
ogy implies cr® y = 2cr^ K = e 2 /h . 30,31 For M 0 S 2 we find 
(J x , y C = 0.71 e 2 /h and since each spin channel contributes 
an equal amount of curvature, this corresponds to 71% 
of the result for the massive Dirac Hamiltonian. 

As discussed in the introduction, measuring the valley 
Hall effect requires generating a carrier imbalance be¬ 
tween the two valleys. Since experimental realizations 
usually involves a gate, we include an overall Fermi level 
ep such that the Fermi level in one valley is £p and in the 
other valley it is £p 5£ due to the presence of polarized 
light. If we assume left-handed polarized light, we can 
then calculate the photoinduced AHC as a function of 5s 
in the K valley, which will depend on the strength of the 
optical perturbation in an experimental setup. Combin¬ 
ing Eqs. (6) and (10) we obtain for sp > A, 


HD 

We can also express this quantity in terms of carrier im¬ 
balance 5n c between the conduction-band edges in the 
two valleys. For a single valley 


7 xk 2 F s 2 f — A 2 
(27t) 2 47 TV 2 h 2 


so that 

2 spSs 5s 2 
4S¥ 

For S£ <C sp we then get 


c 0 C ^ r 

S(Txy ~1l^ S 


-5n r , 


( 12 ) 


(13) 


(14) 


Thus, in this limit the photoinduced AHC becomes linear 
in both the small energy shift 5s and in the small carrier 
imbalance 5n c . An expression similar to this one was 
derived in Ref. 10, but with the Fermi level in the gap. 
That expression can be obtained by setting sp = A in 
Eq. (14). As will be shown in Sec. HIE, the position of 
the Fermi level is minuscule for the intrinsic as well as the 
sidejump contributions, but can have a large effect on the 
skew scattering contribution, which becomes dominant 
for clean samples. 


2. First-principles calculations 

In Fig. 3 we show the intrinsic photoinduced AHC 5cr® y 
from a single spin channel as a function of the single¬ 
valley energy shift fc, for Sp = 0. It should be noted 
that due to the spin-orbit split valence bands in M 0 S 2 , it 
is possible to selectively excite a single spin channel by us¬ 
ing an optical frequency tuned to the transition between 
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FIG. 3. (Color online). Intrinsic photoinduced anomalous 
Hall conductivity of a single layer of M 0 S 2 without gating 
(sf — 0), plotted as a function of the Fermi-level shift Se in 
the K valley. The thin green line corresponds to Eq. (14) 
for the massive Dirac model with the parameters given in 
the text, and the blue dotted line was obtained by evaluating 
Eqs. (5) and (6) from first-principles. 


the topmost valence band and the conduction band. 4,10 
The ab initio calculations are seen to agree very well with 
the model results at the bottom of the valley. When the 
energy shift approaches ~ 0.1 eV the secondary valley 
located between T and K starts to contribute, and the 
model results become unreliable. 

The measurements of the valley Hall effect reported 
in Ref. 10 involved photoexcitation of states in a single 
valley. In that work, the carrier density was estimated 
from photoconductivity measurements, and the photoin¬ 
duced AHC was displayed as a function of carrier density. 
The typical density of photoexcited carriers reported in 
Ref. 10 is of the order of 10 -11 cm -2 , which corresponds 
to an energy shift Se ^ 1 meV; this is far below the point 
where the linear model (14) breaks down. 


III. DISORDERED MOS 2 

Experimental as well as theoretical studies has recently 
demonstrated that sulfur vacancies are the dominant 
source of disorder in M 0 S 2 . 13-18 In the following we will 
thus exclusively focus on the sulfur vacancies and calcu¬ 
late how they affect the VHC. The results will be ana¬ 
lyzed in terms of an impurity averaged unfolded Berry 
curvature, which will be defined below. However, before 
we delve into the ab initio calculations we will briefly 
review the theoretical results for the AHC in a massive 
Dirac model with impurity scattering. 


A. Massive Dirac model 

To obtain a model result for the VHC of disordered 
M 0 S 2 , we use the extrinsic contribution to the AHC of 
the massive Dirac Hamiltonian (7), which has been cal¬ 
culated in the limit of weak and dilute scattering. 12 The 
impurity potential for a given configuration is V(r) = 
EiV^r-Ri), where Ri are impurity sites. The result 
for the intrinsic (0), side jump (SJ) and skew scattering 
(SS) contributions to the VHC of the K valley is then 


e 

2 A 

(15) 

tJ 2h^/(vhkF) 2 + A 2 ’ 

e 2 A 

4 (vhkp) 2 

(16) 

2 h^/(vhk F ) 2 + A 2 

4A 2 + ( vhkp ) 2 


3 (vhk F ) 4 

+ (4A 2 + (vhkp) 2 ) 2 

ss,k e2 (^ 3 )cA ( vhkp ) 4 

xy hx(V 2 ) 2 c (4A 2 + (vhkF) 2 ) 2 ’ 1 ; 

Equation (15) is just Eq. (10) recast in a different form. 
As mentioned in the Introduction, scattering contribu¬ 
tions which are independent of the impurity concentra¬ 
tion x are classified as side-jump, and those which scale 
inversely with x are classified as skew-scattering. In par¬ 
ticular, the second term in Eq. (16) originates from a 
fourth-order expansion of the scattering rate, and could 
thus be regarded as a skew-scattering contribution. Fur¬ 
thermore, the side-jump contribution contains contribu¬ 
tions that cannot directly by ascribed to a coordinate 
shift. 12 

In Fig. 4 we show the VHC as a function of impurity 
concentration at different Fermi level shifts calculated 
within the model with A = 0.86 eV and (V?) c /(V 2 )l = 
—0.4 eV -1 . The VHC converges to the side jump at large 
impurity concentrations, but one should keep in mind 
that the model results are derived under the assumption 
of dilute disorder and weak scattering. The side jump 
correction to the intrinsic part becomes larger for higher 
density of states at the Fermi level and is neglectable 
at very low carrier concentrations. Furthermore, at low 
carrier concentrations, the skew scattering only becomes 
significant at very low impurity concentrations. It thus 
appears that for low carrier concentrations, the intrinsic 
contribution gives a good account of the AHC - even at 
rather low impurity concentrations. 

B. Unfolded band structure and Berry curvature 

In pristine systems, the Berry curvature provides a use¬ 
ful k- resolved measure of the Hall conductivity. For ex¬ 
ample, from Fig. 2 it is clear that states near the valleys 
at K and K r have the largest potential for contributing to 
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FIG. 4. (Color online). Valley Hall conductivity of the dis¬ 
ordered massive Dirac model, Eqs. (15)-(17), evaluated as a 
function of impurity concentration x for three different gating 
levels, with (V?) c /(V?) 2 C — —0.4 eV and A = 0.86 eV. The 
squares at x — 0 correspond to the intrinsic contribution of 
Eq. 15. 


FIG. 5. (Color online). Unfolded spectral function (top) and 
Berry curvature (bottom) of a 3 X 3 M 0 S 2 supercell with a 
single sulfur vacancy. The Fermi level £f = 0 is in the gap 
and the Berry curvature is largely unaffected by the presence 
of the impurity (the Berry curvature for pristine M 0 S 2 , taken 
from Fig. 2, is drawn as a dashed green line). 


the VHC. However, if we would like to know how a given 
distribution of impurities affects the VHC, this picture 
immediately breaks down since the pristine Brillouin zone 
is no longer relevant due to broken translational symme¬ 
try. On the other hand, if a given impurity distribution 
is represented in a supercell, the VHC will still be given 
as a &-space integral of the Berry curvature, but now 
the domain will be the Brillouin zone corresponding to 
the supercell (SBZ), which is not directly comparable to 
the normal Brillouin zone (NBZ). Nevertheless, we can 
expand the supercell states in terms of states in the pris¬ 
tine system and thus unfold the supercell curvature to 
the pristine Brillouin zone. For a general band quantity 
cln( K) defined in SBZ we can thus define the unfolded 
quantity a^ (k) in the NBZ as 

4“)(k) = Y, |(VK|nk)| 2 ajv(K), (18) 

N 

where K is the SBZ crystal momentum, which is related 
to k by translation of a supercell reciprocal lattice vector. 
Brillouin zone integrals can be written in terms of the 
unfolded function since 


The method of unfolding k -space quantities has pre¬ 
viously been applied to band structures of disordered 
systems 27,28 and more recently to unfolding the Berry 
curvature. 29 In the case of band structures the object of 
interest is the spectral function 

A w (u,k) =^KVK|nk)| 2 «5( W - ej vK). (20) 

nN 

The Berry curvature is somewhat more complicated 
since, a naive application of Eq. (18) leading to 

fi (u) (k) = ]T \(NK\nk)\ 2 f NK Q N (K), (21) 

nN 

becomes gauge dependent. In Ref. 29 Eq. (21) it was 
shown h how to generalize Eq. (18) to obtain an ex¬ 
plicitly gauge invariant quantity. In the present work we 
have applied the gauge invariant expression for all calcu¬ 
lations. 


C. Sampling impurity configurations 


A= [ dKVay(K) 

J SBZ N 

= [ dK [ dkV|(VK|nk)| 2 aAr(K) 

J SBZ J NBZ nN 

= [ dk EXV k )> ( 19 ) 

J NBZ 


where a complete set of NBZ states were inserted in the 
second line and it was used that (iVK|nk) is only non¬ 
vanishing if k downfolds to K. 


In Fig. 5 we show the unfolded spectral function and 
Berry curvature of a periodic structure obtained as 3 x 3 
M 0 S 2 unit cell with a single sulfur vacancy. The vacancy 
is seen to introduce both occupied and unoccupied states 
in the gap, but the Berry curvature is largely unaffected 
by such rather localized states. This system is a example 
of a x = 1/18 = 0.56 impurity concentration, but is not 
necessarily representative of this disorder concentration 
in general. As it turns out the Berry curvature is largely 
insensitive to the impurity configuration as long as the 
Fermi level is in the gap. However, the situation changes 
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FIG. 6. (Color online). Same as Fig. 5, except that the Fermi 
level has been shifted to the conduction band, where the im¬ 
purities induce large modifications of the Berry curvature. 

dramatically when the Fermi level is shifted to the con¬ 
duction band. This situation is shown in Fig. 6. The 
blurred features in the unfolded spectral function is as¬ 
sociated with scattering states and are accompanied by 
spiky feature in the Berry curvatures. From the spec¬ 
tral representation of the Berry curvature Eq. (Al), it 
is clear that such features arise whenever occupied and 
unoccupied states come close to the Fermi level. 

In experiments with M 0 S 2 transistors the Fermi level is 
typically controlled with a gate voltage. 10 Furthermore, 
exfoliated M 0 S 2 usually exhibits an intrinsic n-doping, 
which is attributed to Re impurities. 19 For these reasons, 
we have chosen to pin the Fermi level to the conduction 
band in the following. This will facilitate the comparison 
with experiments as well as model calculations. 

In order to perform a faithful calculation of the conduc¬ 
tivity at a given impurity concentration, one then has to 
average over a large number of impurity configurations. 
Such a procedure requires a high number of simulations in 
large unit cells and is not feasible with standard ab initio 
methods. To proceed we construct an effective Hamilto¬ 
nian based on ab initio DFT calculations and Wannier 
functions. Given an impurity concentration, we thus ran¬ 
domly generate an impurity configuration and construct 
the corresponding Hamiltonian as described in Appendix 
B. For a given vacancy concentration, we perform cal¬ 
culations for ^ 1000 systems with randomly drawn im¬ 
purity configurations in a 6 x 12 supercell. The smallest 
impurity concentration considered was x = 1/216, where 
we needed a larger supercell of 12 x 18. As a first check 
of the method, we calculate the optical conductivity of 
M 0 S 2 at various S vacancy concentrations. The results 
are shown in Fig. 7 and as expected the impurities sim¬ 
ply introduce a broadening of the spectra. We should 
note that the Wannier functions used to construct the 
tight-binding Hamiltonian was obtained by disentangling 
bands up to 3.0 eV above the conduction band minimum 
and the present calculation is thus only expected to be 
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FIG. 7. (Color online). Configuration-averaged optical ab¬ 
sorption spectrum of M 0 S 2 at various impurity concentration. 


reliable up to 3.0 eV above the pristine absorption edge. 


D. Valley Hall conductivity 

The lowest impurity concentration we have considered 
is x = 0.005. Here the configurational averaged unfolded 
Berry curvature becomes very spiky and is not particu¬ 
larly informative. In Fig. 8 we show the configurational 
averaged unfolded Berry curvature for 4 different inter¬ 
mediate impurity configurations for £p — A = 0.24 and 
£p — A = 0.24. In general we observe that the peaks in 
the Berry curvature becomes enhanced and broadened, 
while the curvature retains its qualitative features. This 
tendency is maintained over a range of impurity concen¬ 
trations (x = 0.02 - x = 0.1) and gives rise to a impurity 
independent increase in VHC. We will identify this with 
the side jump corrected VHC and the Berry curvatures in 
Fig. 8 thus comprises a measure of the k -space resolved 
side jump scattering. From the unfolded spectral func¬ 
tion it is clear that the sulfur vacancies has the effect of 
lowering the overall potential, such that the conduction 
bands are lowered with respect to the fixed Fermi level 
at high impurity concentrations. However, for interme¬ 
diate dopings the exact position of the Fermi level does 
not have a large effect on the VHC. In fact, raising the 
Fermi level with respect to the conduction band tends to 
lower the VHC, since more curvature from the conduc¬ 
tion band will be included and this has the opposite sign 
of the dominating contribution from the valence bands. 
At larger impurity concentrations (x = 0.125) the Berry 
curvature becomes more ” smeared” and will eventually 
average to zero. At this point the system is so strongly 
perturbed that it cannot be analyzed in terms of scatter¬ 
ing events. 

In Fig. 9 we show the VHC as a function of impurity 
concentration at different Fermi level shifts. At low im- 























FIG. 8. (Color online). Impurity averaged unfolded spectral function and Berry curvature of disordered M 0 S 2 with the Fermi 
level fixed in the conduction band. The dashed green lines shows the result for pristine M 0 S 2 . 


purity concentrations we observe a divergent behavior of 
the VHC. Comparing with Eqs. (15)-(16) we identify this 
with skew scattering processes. In terms of the unfolded 
curvature, this is the ” spiky” regime, where the averaged 
curvature largely looses the qualitative features of the 
pristine system. We note that as the impurity concen¬ 
tration x is decreased, it becomes progressively harder to 
converge the results. This is due to the fact that the stan¬ 
dard deviation for a given supercell size increases when 
x becomes small, while at the same time we need very 
large supercells in order to perform calculations for small 
x. 

At intermediate impurity concentrations the VHC 


reaches a flattened plateau, which we identify as the side 
jump regime. In this regime skew scattering is insignif¬ 
icant and the pristine result receives a small correction 
which is nearly independent of impurity concentration. 
Comparing with Fig. 4, we see that there is good qual¬ 
itative agreement between the model and the ab initio 
results. A major difference is the fact that the ab initio 
VHC decreases at high impurity concentrations, whereas 
the model converges towards the side jump result. The 
reason is that the model results were derived under the 
assumption of dilute impurity concentration, and cannot 
be applied in this regime. Moreover, at large Fermi level 
shifts the deviations from the Dirac model begins to be 
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FIG. 9. (Color online). Calculated VHC as a function of im¬ 
purity concentration x. We show the result for three different 
positions of the Fermi level with respect to the conduction 
band minimum e c . 
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FIG. 10. (Color online). Photoinduced AHC as a function 
of Fermi level shift 5s between the valleys. States in the K 
valley are assumed excited corresponding to application of a 
left-handed polarized field. The dotted lines are the intrinsic 
5cr the dashed lines are 5a^ y + 5a^ y , and the solid lines 
are 5a^, y + 5a^ y + 5a^ y . Note the large dependence on Fermi 
level when skew-scattering is included. 


important. In the present system the most important 
effect is the contribution from the secondary conduction 
band minimum between T and K. 

When the Fermi level comes close to the conduction 
band minimum, two crucial features can be observed 
from Figs. 4 and 9. First, the side jump correction 
becomes small such that the intrinsic VHC becomes a 
good estimate at intermediate impurity concentrations. 
Second, the critical impurity concentrations, where skew 
scattering starts to dominate, becomes rather small. For 
Sp — = 0-04 this concentration is x ~ 0.005 and will 

be even smaller when the Fermi level moves closer to the 
band edge. 


except that the sign has been changed. This was al¬ 
ready noted in Ref. 10 For the VHC, the side jump only 
contributes a minor part to the scattering independent 
VHC. However, for the derivative, which is the relevant 
quantity for the photoinduced AHC, the side jump plays 
a dominant role. At small impurity concentrations, the 
skew scattering contribution will dominate. The skew 
scattering part of the VHC scales as 1/x and so does 
its contribution to the photoinduced AHC. Interestingly, 
the skew scattering contribution also has significant de¬ 
pendence on the position of the Fermi level and could 
explain the large dependence on gate voltage observed 
experimentally. 10 


E. Photoinduced anomalous Hall conductivity 

As mentioned previously, the experimentally measured 
quantity is the photoinduced AHC, which we express as 
the sum of the two VHCs. For small carrier imbalance 
this will be Sa xy (sF,Ss) = davH{^)/ds\ £=£F 5s, where 5s 
is the shift in Fermi level corresponding to the charge 
imbalance. A first principles evaluation of this expres¬ 
sion including impurities is made difficult by the fluc¬ 
tuations involved in the configurational averaging pro¬ 
cedure. However, the good agreement between the ab 
initio calculations of the VHC (Fig. 9) with the mas¬ 
sive Dirac model (Fig. 4), suggests that we can obtain 
a reliable estimate of the photoinduced AHC from the 
Dirac model. For the skew scattering contribution we 
take (V?) c /(V?)c = —0.4 eV as in Fig. 4. Note that this 
value, and in particular the sign, is a non-trivial ab initio 
result obtained from the VHC calculations. 

In Fig. 10 we show the photoinduced AHC calculated 
from Eqs. (15)-(17) with x = 0.01. First of all it should 
be noted that the inclusion of side jump scattering gives 
a result, which is very close to the intrinsic contribution, 


IV. CONCLUSIONS 

A priori it is not clear that the intrinsic VHC provides 
a good descriptor for the VHC in a M 0 S 2 transistor setup, 
since the Hall conductivity is expected to diverge in the 
clean limit as a consequence of skew scattering. 10 Here 
we have used first principles calculations to investigate 
the effect of sulfur vacancies at different impurity con¬ 
centrations. The influence of disorder was analyzed in 
&-space in terms of the unfolded Berry curvature and we 
have shown that the side jump regime appears as a con¬ 
centration independent enhancement of the Berry curva¬ 
ture. The skew scattering introduces divergences in the 
Berry curvature and the unfolded Berry curvature be¬ 
comes spiky and irrelevant. Nevertheless, we were able 
to converge the VHC calculations in the skew scatter¬ 
ing regime, and recover the expected divergent behav¬ 
ior. The ab initio calculations show qualitative agree¬ 
ment with model calculations based on a massive Dirac 
Hamiltonian and allow us to extract a non-trivial value 
of the skew scattering potential (V?) c /(V?)^ = —0.4 eV. 
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The calculations allow us to estimate a critical impurity 
concentration x where skew scattering starts to domi¬ 
nate. For £f = 40 meV (with respect to the conduction 
band minimum) we find that x ~ 0.005 and below this 
point the intrinsic VHC becomes a poor descriptor. 

A comparison with experiments 10 indicates that in¬ 
deed, the intrinsic VHC cannot be applied as a descriptor 
for the photoinduced valley Hall conductivity. As previ¬ 
ously noted, the side jump contribution changes the sign 
of photoinduced AHC and we have shown that the skew 
scattering contribution is a likely explanation for the 
large gate dependence observed experimentally. 10 How¬ 
ever, a reliable estimate of the skew scattering contribu¬ 
tion requires knowledge of the impurity concentrations 
in the samples investigated, which is not presently avail¬ 
able. It would be very interesting to perform measure¬ 
ments of the photoinduced AHC on M 0 S 2 samples with 
different impurity concentrations in order to unravel the 
roles played by side jump and skew scattering. 

We have implicitly assumed a low temperature regime 
and therefore not discussed the effect of phonons. The 
effect of phonons on the longitudinal mobility in M 0 S 2 
has been analyzed thoroughly, 32 but the influence on the 
transverse conductivity has so far not been considered. 
Furthermore, monolayer M 0 S 2 have been shown to ex¬ 
hibit strong excitonic effects 33,34 due to poor screening 
in two-dimensional materials. The charge imbalance uti¬ 
lized in the experimental realization originates from opti¬ 
cally generated electron-hole pairs, and if the Fermi level 
is close to the conduction band edge these effects could 
potentially severely limit the carrier mobility. We will 
leave these issues for future studies. 
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Appendix A: Calculational details 

The calculations in the present work was performed 
with the tight binding method using parameters obtained 
from ab initio density functional calculations and pro¬ 
jected Wannier functions. 


1. Ab initio calculations and construction of 
Wannier orbitals 

The ab initio density functional theory calculations 
were performed with the pwscf code from the Quantum 
Espresso package, 35 using the PBE functional. Norm- 
conserving pseudopotentials were used, and the calcula¬ 


tions were carried out with a planewave cutoff of 100 Ry. 
The lattice parameter of M 0 S 2 was set to the experi¬ 
mental lattice constant of 3.16 A, and 12 A was used to 
separate the periodically-repeated images. All calcula¬ 
tions were performed in a non-collinear spin framework, 
with fully relativistic pseudopotentials. 

After converging the Kohn-Sham electronic structure, 
the valence and low-lying conduction Bloch bands were 
converted into projected Wannier functions using the 
Wannier90 code package. 36 For M 0 S 2 the projected Wan¬ 
nier orbitals were constructed using sulfur p states and 
Mo d states. The sulfur s states were included in the 
ab initio calculations, but the low lying 8-like Bloch 
bands were excluded from the wannierization. Unoccu¬ 
pied states were included by disentangling 37 bands up 
to 3.0 eV above the conduction band minimum. Fi¬ 
nally the Wannier functions were used to construct the 
Kohn-Sham Hamiltonian in a tight-binding basis Hjj( R), 
where i, j denotes orbital indices within the unit cell and 
R is a lattice vector. The set of lattice vectors included 
were defined by the Wigner-Seitz cell corresponding to 
the applied ab initio k -point mesh. For example, in pris¬ 
tine M 0 S 2 with a 8 x 8 k -point mesh, we have 22 orbitals 
(Mo d and S p) and 64 lattice vectors (some of which are 
equivalent). 

2. Tight-binding calculations 

The majority of calculations in present work are tight- 
binding calculations with parameters obtained from a 
Wannier representation of the Kohn-Sham Hamiltonian 
Hij (R) . In a tight-binding framework, the calculation 
of band structures from Hjj (R) is of course equivalent 
to the standard Wannier interpolation. 37 At the sam¬ 
pled set of ^-points, Hjj (R) will thus yield the calculated 
Kohn-Sham eigenvalues and between the sampled points 
it smoothly interpolates. Similarly, a rigorous Wannier 
interpolation scheme can be constructed for the AHC, 38 
but this quantity cannot be calculated exactly in a bare 
tight-binding framework since the information contained 
in Hjj (R) is not enough to evaluate the Berry curvature. 

To see this explicitly we will briefly state the relevant 
expressions below. The starting point is the Berry curva¬ 
ture in its spectral representation where it can be written 

^a^(k) = ^(/nk - /mk) (Al) 

m,n 

^ (^mk| V a H (k) ju-nk) {u n k| V pH (k) ju-mk) 

(^nk ^mk)^ 

with V a = d/dk a . We let cfj denote at set of localized 
orbitals and expand the Bloch states as 

iV’nk) = J^C'mklXik) = E^ k ^ R+t *)|^R) (A2) 

i iR 

where t j = (4>io\^\4>io)• Note that the inclusion of t j 
is purely a matter of convention. However, the present 
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convention will prove highly convenient below. Matrix and in terms of these, the matrix elements appearing in 

elements of the Bloch Hamiltonian and their gradients Eq. (Al) becomes 

can now be written as 

Hi, k = (Xik|-ff |Xjk) = E e* k ' ( R-ti + V Hij (R), (A3) 

R 

V k if ijk = i R - t< + (A4) 

R 


Jamn (k) — (^mk |rt n k) 

= E^kCni k e ik - (R - ti+t ^(^o|e ik -W QJ ff(k)e- k -% R ) 

ijR 

= -*E ^kC'nike ik -( R - t ‘ +t ^ <^o|[r, tf] |^r) 
t?R 

= M E G mik^nike* k ' (R_ti+ti y^io|f|m')(m'l^l^iR) - (v>io\H\ipi R ,){<pi R ,\r\(p jR )] 
r/ZRR' 

= M ^ C , mikC , njke lk '^ R+R [r^R/i^jR - 7U/R'(UjR + R'^RO^j)] 

r/ZRR' 

= H E CmikCnikkiiicHyic - H m r ljk \ + i E C , ^ ik C , nj - k Re ik -( R - t » +t ^fl iiR 

ijl ij R 

= ^ ^ ^raik^njk V a -H^k + ^ ^ ^ ^raik^C/k r i^k)-f^Zjk T/^Zk(U^Zj r Zjfk)] * (A5) 


U ijl 


In the present work we have made the diagonal ap¬ 
proximation where r= ^Ro^ijU and we thus use 
4mn(k) = J2ij C* l(k Cn jk V a i?,; jk . With this approxi- 
mation, the problem can be mapped exactly to a tight- 
binding calculation with parameters obtained from the 
Kohn-Sham Hamiltonian in a basis of Wannier functions. 
This turns out to be an excellent approximation. The 
Berry curvature calculated with this method cannot be 
distinguished with the naked eye from ab initio calcu¬ 
lations (obtained by Wannier interpolation) and anoma¬ 
lous Hall conductivities calculated within the diagonal 
approximation differ from ab initio results by less than 
2 %. 


3. Mapping the supercell onto the normal cell 

Unfolding band structures and curvatures involve the 
calculation of matrix elements between pristine and su¬ 
percell systems. A localized basis set allow us to per¬ 
form the unfolding without direct reference to the pris¬ 
tine system. 28 

We denote pristine band, orbital indices, and crystal 
momentum by n, i, k respectively and super cell band, or¬ 
bital, and crystal momentum indices by TV, /, K respec¬ 
tively. For the present purpose we will use r for pristine 
lattice vectors and R for supercell lattice vectors. We 


can thus consider the matrix element 

(V’nkIV’iVK) = E < ^nikCw/Ke _ * k ' (r+ti)+ * K ' (R+t/) (^ir|^/ R ) 

il rR 

(A6) 

Following Ku et al. 28 we introduce a map that uniquely 
identifies orbitals in the supercell with corresponding 
orbitals in the pristine system. The map thus takes 
I —>• r'(/),z'(/) and we have 

(</Mc|^/r) = (^ir|^z'(/)R+r / (7)) = fe(7)Wr'(/)- (A7) 

The matrix element can then be written as 

(V’nkIV’iVK) = E C '^(/)kCiV/Ke- lk - [r ' (/)+R+t '(^) 1+iK - [R+ ^ 

in 

= E C , ;i'(/)kC'i V 7Ke- k -t r '( / ) +t ^^) +<K -t^ K[k] , 

(A8) 

where [k] is the set of crystal momenta that downfolds 
to K. 

The idea of a map allows one to only work with the su¬ 
percell system and avoid the explicit calculation of over¬ 
lap matrices. However, the procedure does require that 
the supercell system considered is somewhat similar to 
the pristine reference system and becomes ill-defined if 



12 




FIG. 11. (Color online). Unfolded band structure and Berry curvature of M 0 S 2 in a 3 X 3 supercell with two sulfur vacancies 
at next nearest neighbor sites. Left: ab initio calculation. Right: calculation from effective Hamiltonian constructed from 
a pristine 3x3 supercell with impurities added in the form of an influence Hamiltonian obtained from a single impurity 
calculation. The dashed lines show the Berry curvature of pristine M 0 S 2 with the Fermi level at the indicated position. 


there is not a unique way of relating orbitals in the two 
systems. Furthermore, it is important to construct the 
tight-binding Hamiltonian from projected Wannier func¬ 
tions as opposed to maximally localized Wannier func¬ 
tions, since the latter may differ significantly in otherwise 
similar systems. 

In Figs. 5 and 6 we showed examples of the band 
structures and Berry curvature of M 0 S 2 in a 3 x 3 unit 
cell with a single sulfur vacancy, unfolded onto the normal 
BZ. 


Appendix B: Effective Hamiltonians 

Here we briefly summarize the construction of effective 
Hamiltonians as proposed in Ref. 27. In a tight binding 
framework, the effective Hamiltonian with N impurities 
is constructed as 

HfJ f (R) = -ffff(R) (Bl) 

N 

+ Y, P{I, J , r a , R)ff', (/) ;, (J) (!•'(/) - r„ r'( J) + R - r s ) 

S = 1 

where /, J, R denotes a supercell orbitals and lattice vec¬ 
tors respectively, r is a normal cell lattice vector and 
i',r' are maps from the supercell orbitals and lattice to 
the normal cell. r s denotes the position of impurity s. 
The influence Hamiltonian is given by 

ri, r 2 ) = r 2 ) - ff?.(r 2 - n), (B2) 

where ri,r 2 ) is constructed with a map from the 

impurity to the normal cell. Since this is expressed in 
a basis of normal cell lattice vectors H Imp and there¬ 
fore H' is not periodic in simultaneous translations of iq 
and T 2 . We have used the partition function of Liu and 


Vanderbilt 39 

P(I, J,r„R) = e~ Wd ^ 8 (B3) 

d = (|r'(J) - r s | + |r'(J) + R - r s \)/2. 

with d c = 9.0 A. 

For applications of the method we have performed ab 
initio calculations of the pristine and impurity systems 
and constructed H irnp and H° using projected Wannier 
functions. H sc is then constructed as a straightforward 
repetition of H°. 

As a non-trivial test of the method we have performed 
an ab initio calculation of a 3 x 3 unit cell of M 0 S 2 
and with sulfur vacancies at two next-nearest neighbor 
sites. We have then constructed the same system with 
from the effective Hamiltonian: first we construct obtain 
the tight binding model of a 3 x 3 unit cell of M 0 S 2 by 
repeating the tight-binding Hamiltonian obtained from 
a calculations of pristine M 0 S 2 (normal unit cell). We 
have then constructed the influence Hamiltonian (B2) 
from the pristine calculation and a calculation of M 0 S 2 
5 in a 3 x 3 unit cell with a single sulfur vacancy. The 
influence Hamiltonian is then added to the 3x3 tight- 
binding Hamiltonian at the two nearest neighbor sites to 
obtain the effective Hamiltonian of a system with two 
impurities. The unfolded bands and curvature of ab ini¬ 
tio and effective Hamiltonian calculations are shown in 
Fig. 11. The unfolded spectral function is nearly in¬ 
distinguishable in the two case. However, the unfolded 
Berry curvature is highly sensitive to the exact positions 
of bands near avoided crossings and therefore exhibits 
larger deviations in the two methods. Nevertheless, the 
Berry curvature obtained from the effective Hamiltonian 
reproduces the main qualitative features (for example a 
vanishing contribution in the vicinity of K and —K) and 
we believe that the configurational averaged Berry cur¬ 
vature obtained from the effective Hamiltonian methods 
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provides a reliable measure of the effects of disorder on the curvature and VHC. 
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